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lfa  M#yi»tf<^  mt«  q|  tiMi  >rt  lo  co^ting  t»cft»>»loty 
«g4  ;w»m-:a  Sithods  brings  increasingly  difficult 
#I^*HJIU  («oM«u  iBW«tj»tlo*.  Som  of 

thays  pcoblams,  -Ip.  turn,  challenge  tha  capabilities  of 
Avtiifcibg  ona«rio«l  ttdtaipM  iad  stinulate  develop¬ 
ment:  of  new  oaas'la  1981  the  Office  of  Msvel  Research 
pet  forth  0  call  to  two  national  laboratories  for  for¬ 
th#*  lwolwint  in  a  longstanding  problse  In  nave 
propagation  theory:  tha  behavior  of  a  nonlinear  pulse 
incident  ou  a  caustic  surface.  This  had  been  an  unre- 
solved  issue  in  tha  aircraft-generated  sonic  boon 
studies  of  the  1970s •  ft  is  currently  an  issue  in  tbs 
oceanic  aselsnr  standoff  problea.  Experimentation  with 
existing  methods  for  nonlinear  hyperbolic  system  sad 
shock  waves  prompted  development  of  the  numerical 
technique  preseaced  here,  lech  la  its  present  form  and 
In  future  reflnenents  this  technique  can  be  of  value 
in  other  .  areas,  tramples  of  naval  interest  are  Inroad* 
band  pulse  propagation,  surface  dad  internal  eaves, 
end  ocean  circulation. 


Executive  Summary 


,>  We  present  first-  and  second-order  upwind  schemes  em¬ 
ploying  a  numerically  calculated  characteristic  speed 
direction  and  combine  them  into  a  simple  monotonicity 
preserving  hybrid  scheme  using  the  method  of  flux  cor¬ 
rection.  The  first-order  scheme  is  constructed  to 
maintain  accuracy  at  flow  reversal  points.  The  hybrid 
scheme  computes  a  provisional  update  from  the  first- 
order  scheme,  and  then  filters  the  second-order  cor¬ 
rections  to  prevent  occurrence  of  new  extrema.  We  de¬ 
rive  analytic  solutions  for  a  developing  N-wave  shock, 
and  compare  computed  versus  analytic  results  for  two 
different  N  waves  and  for  a  third  case  involving 
linear  advection  of  a  square  wave.  Results  are  given 
with  and  without  the  second-order  correction.  The 
second-order  results  are  always  superior  to  first- 
order  results,  with  the  most  dramatic  difference  oc¬ 
curring  for  the  case  of  linear  advection.  The  results 
suggest  that  higher  order  differences  could  be  substi¬ 
tuted  in  the  hybrid  scheme  to  reduce  truncation  error 
even  further 
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Illustrations 

Figure  l.  Symmetric  N-wave  shock  development.  11 

Solid  line:  analytic  solution;  dotted 
line:  numerical  solution.  Dot-dash 
line:  domain  of  initial  condition. 
First-order  results  are  given  in  b-d; 
second-order  in  e-g. 

Figure  2.  Drifting  N-wave  demonstrating  stability  12 
against  expansion  shocks.  The  character¬ 
istic  speed  vanishes  where  the  amplitude 
is  35  percent  of  initial  maximum.  For 
legend  see  Figure  1.  The  step-like 
structure  at  step  601  is  a  result  of 
periodic  boundary  conditions. 

Figure  3.  Linear  advection  of  a  square  wave.  Solid  13 
line:  exact  solution;  dashed  line: 
numerical  solution.  First-order  results 
appear  in  b-d;  second  in  e-g. 


A  Second  Order  Upwind  Flux  Method 
for  Hyperbolic  Conservation  Laws 


1.  Introduction 

The  development  of  flow  discontinuities 
in  hyperbolic  systems  presents  substan¬ 
tial  difficulties  to  finite  difference 
solution  algorithms.  First-order  upwind 
differences  taken  in  the  direction  of 
the  characteristic  speed  alleviate  art¬ 
ificial  oscillations  near  a  shock  at 
the  expense  of  adding,  implicitly  or 
explicitly,  artificial  dissipation.  The 
flux  correction  approach  of  Zalesak 
(1979)  addresses  this  problem  through 
construction  of  a  hybrid  scheme  that 
results  in  high-order  centered  differ¬ 
ences  being  used  where  the  flow  is 
smooth,  and  one-sided  differences  being 
used  near  discontinuities. 

Although  this  method  is  capable  of  pro¬ 
ducing  impressive  results  for  initial 
conditions  containing  shocks  (Zalesak, 
1981),  it  generates  shocks  too  rapidly 
in  some  cases  where  they  evolve  from 
continuous  initial  conditions.  Here  we 
demonstrate  that  the  flux  correction 
technique  enables  the  use  of  spatially 
high-order  upwind  differencing  while 
maintaining  proper  monotonicity  con¬ 
straints.  We  show  that  a  hybrid  scheme 
consisting  of  first-  and  second-order 
upwind  differences  in  tandem  results  in 
excellent  agreement  between  analytic 
and  numerical  solutions  for  a  steepen¬ 
ing  N-wave  shock — a  familiar  phenomenon 
in  nonlinear  acoustics.  The  scheme  pre¬ 
sented  here  is  under  continued  investi¬ 
gation  to  remove  expansion  shocks 
(time-reversed  shocks  forbidden  by  the 
second  law  of  thermodynamics)  and  a 
mild  amplification  of  long  wavelength 
modes . 

We  wish  to  distinguish  between  the  ap¬ 
proach  taken  here  and  that  taken  in 
contemporary  versions  of  the  scheme  of 


Godunov  (1959).  These  schemes  are 
reviewed  by  van  Leer  (1981)  and  by 
Harten,  et  al.  (1983).  These  schemes 
view  flow  variables  as  piecewise  con¬ 
stant,  so  that  a  continuous  variation 
becomes  the  limit  of  vanishingly  small 
stepwise  variations.  Time  advancement 
is  achieved  by  invoking  exact  or  ap¬ 
proximate  solutions  of  the  individual 
noninteracting  shock  problems  (i.e., 
the  Riemann  problem  for  Burgers'  equa¬ 
tion)  for  a  small  increment  of  time  St. 
A  new  stepwise  representation  is  then 
recovered  by  spatial  averaging.  This 
approach  has  advantages  when  the  flow 
consists  of  clean,  well-formed  shocks. 
However,  in  regions  where  the  flow  is 
dominated  by  linear  advection,  the  re¬ 
sults  are  highly  dissipative,  as  with 
the  first-order  upwind  method. 

Here  we  present  the  first-  and  second- 
order  upwind  methods,  derive  an 
analytic  model  for  a  steepening  N-wave 
shock,  and  compare  numerical  versus 
analytic  results  for  this  model  and  for 
the  linear  advection  of  a  square  wave. 
The  first-order  method  differs  from 
others  in  current  use  by  employing  a 
characteristic  speed  calculated  in¬ 
ternally  from  point  values  of  flux  and 
the  variable  of  integration  (density  or 
velocity) .  It  is  also  independent  of 
any  specific  model  such  as  Burgers' 
equation . 

In  Section  2  we  present  equations  and 
conservation  laws  that  motivate  our 
approach.  Section  3  derives  analytic 
models  for  validation  of  the  numerical 
results.  Section  4  presents  the  numeri¬ 
cal  methods  to  be  used,  and  shows  the 
first-order  scheme  to  be  monotone  for 
all  flow  conditions  as  long  as  the 
Courant  condition  is  satisfied.  Section 
5  then  compares  analytic  and  numerical 


1 


results  for  the  second-order  corrected 
scheme  and  for  the  first-order  scheme 
alone . 


~ y*  pmdx  =  0,  m  *  1,  2,  3  ... 


(4) 


2.0  Equations  and  Conservation  Laws 

We  are  concerned  with  hyperbolic  sys¬ 
tems  which  can  be  modeled  by  the  scalar 
initial  value  problem 

P*  +  <J>(p)„  =  0,  (la) 

p(x,0)  =  PQ(x),  (lb) 

where  the  flux  0  is  a  smooth  function 
of  p. 

The  variable  <p  in  various  applications 
may  play  the  role  of  a  mass  density,  a 
density  perturbation  (as  in  the  non¬ 
linear  acoustics  problem  which  motivat¬ 
ed  this  work.),  or  a  velocity  variable. 
Thus,  in  general,  p  may  take  on  posi¬ 
tive  or  negative  values. 

The  method  developed  here  is  intended 
for  systems  more  general  than  Eq.  (1), 
with  source  terms  and  coupling  to  other 
variables.  For  the  moment,  however,  we 
use  well-known  properties  of  Eq .  (1) 

including  analytic  solutions  as  a  test 
of  the  numerics.  Equation  (la)  may  be 
expressed 

Pt  +  u(p)  p  =  0,  (2) 

v  A 

where  u(p)  =  0P  (3) 


is  known  as  the  characteristic  speed. 
Equation  (2)  implies  preservation  of 
raonotonicity :  pt  *  0  at  local  extrema. 

It  is  also  worth  noting  that  prior  to 
shock  formation,  Eq.  (la)  implies  under 
a  variety  of  boundary  conditions  that 
all  moments  of  the  distribution  p  are 
constant : 


Depending  upon  the  form  of  the  flux  and 
the  initial  condition,  Eq.  (2)  may  ad¬ 
mit  analytic  solutions  up  to  the  forma¬ 
tion  of  a  shock.  These  are  found  as 
follows.  Equation  (2)  implies  constancy 
of  p  along  the  trajectory 

Xt  =  u(p).  (5) 


Both  xt  and  u  are  constant  on  this  tra¬ 
jectory,  so  the  value  of  p  at  (x,t)  may 
be  found  by  mapping  back  to  the  initial 
condition : 

p(x,t)  =  P0(x  -  t  u(p)).  (6) 

Analytic  solutions  are  generated  by 

solving  Eq.  (6)  for  p  in  terms  of  x  and 

t.  In  an  example  to  be  given  in  Section 

3,  u(p)  is  taken  to  be  linear  and  p  (x) 

o 

quadratic.  Then  p  (x,t)  emerges  as  the 
solution  of  a  quadratic  equation. 

After  sufficient  time  the  characteris¬ 
tics  (Eq.  (5))  for  neighboring  points 
of  an  initial  distribution  may  cross, 
leading  to  mathematical  solutions  which 
are  multivalued.  These  mathematical 
solutions  still  maintain  monotonicity 
and  constancy  of  all  moments  in  Eq . 
(4),  but  loose  physical  meaning.  Physi¬ 
cally,  an  arbitrarily  small  amount  of 
viscosity  becomes  dominant  as  |PX|  aP~ 

proaches  infinity,  and  a  shock  is  form¬ 
ed.  From  that  time  forward,  only  the 
fundamental  moment  (m=l  in  Eq.  (4))  is 
preserved.  This  corresponds  to  conser¬ 
vation  of  mass  or  momentum,  depending 
upon  the  interpretation  of  the  variable 
p.  To  illustrate  this  point,  we  will 
derive  the  Rankine-Hugoniot  relation  as 
follows.  Consider  a  discontinuous  pro¬ 
file  with  constant  left  and  right 
states : 


p  =  PL  +  (pR  "  PL^  H(x-Vt)  (7) 


4>(p)  =  <)>L  +  (4>r  -  <t>L)  H(x-vt) , 


where 
H(x)  = 


1,  x>0 
0,  x<0  , 


(8) 

(9) 


and  v  is  the  shock  propagation  speed. 
Substituting  Eqs. (7)-(9)  into  Eq .  (la) 
and  equating  coefficients  of  H(x)x  - 

6  (x)  gives  the  Rankine-Hugoniot  relation 


v(»L  -  PR)  ■  *L  '  *R  • 


(10) 


Generalizing  Eqs.  (7)-(10),  a  propaga¬ 
tion  speed  could  be  defined  by  im¬ 
posing  conservation  of  the  quantity  pm 
across  the  shock.  Eqs.  (la)  and  (3) 
give 


/P 

u(p)  pm"  dp  =  0,  (11) 

and  Eqs.  (7)— (10)  generalize  to 

(m  m\ 

PL  -  PR  / 

rpl  , 

=  ml  pm"  u(p)  dp  .  (12) 

PR 


Unless  u  *  constant,  corresponding  to 
linear  advection  at  a  constant  velo¬ 
city,  preservation  of  each  moment  would 
require  a  different  propagation  speed. 
Conservation  of  mass  (or  momentum,  de¬ 
pending  upon  the  interpretation  of  p) 
is  absolutely  essential,  so  Eq.  (10) 


wins  out  in  the  determination  of  the 
shock  speed. 

It  is  important  to  note  that  flux  con¬ 
servation  (Eq.  (la))  is  sufficient  to 
yield  the  proper  shock  speed.  This  will 
be  exploited  in  the  numerical  scheme  to 
be  constructed  below.  In  particular, 
the  numerics  are  cast  in  flux  conserva¬ 
tive  form,  with  no  distinction  between 
shocks  and  steep  but  continuous  pro¬ 
files.  The  emergence  of  a  shock  with 
proper  speed  from  a  steepening  continu¬ 
ous  profile  then  serves  to  validate  the 
method . 

3.0  Analytic  Models 

In  this  section  we  present  two  analytic 
models  to  be  used  as  benchmarks  for  the 
numerical  algorithms.  The  first  is  a 
folded  parabolic  lobe  (approximating 
one  cycle  of  a  sine  wave)  that  evolves 
nonlinearly  to  an  N-wave  shock.  The 
second  is  the  linear  advection  of  a 
square  wave  in  a  constant  velocity 
field — a  problem  posed  by  Boris  and 
Book  (1973)  and  used  frequently  to  dem¬ 
onstrate  the  flux  correction  technique. 

Let  us  consider  an  inviscid  Burgers' 
equation  cast  in  a  moving  frame  of 
reference : 


pt  +  (f  pZ  '  V}x  =  °’  (13) 

where  a  and  u  are  constants.  The  u 
o  o 

term  is  added  to  demonstrate  the  numer¬ 
ics'  robustness  against  formation  of 
nonphysical  "expansion  shocks."  Let  the 
initial  condition  be 


P0  (x)  =  max  (0,  bx  (xq-x)),x»0(14) 
and 

p0  (-x)  =  "  P0(x)-  (15) 

The  result  is  a  pair  of  antisymmetric 
parabolic  lobes  shown  in  Figure  1. 
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Within  each  lobe,  an  analytic  solution 
can  be  found  from  Eq.  (6).  From  Eqs. 
(13)  and  (la). 


drops  discontinuously  to  zero.  The 

shock  location  x  can  be  determined  an- 
s 

alytically  by  setting 


U  =  ap  -  UQ  , 

so  the  equation  to  be  solved  is 

p  =  xq  (x  -  t(ap-uc)) 

-  (x  -  t(ap-uQ))2]  . 


The  solution  which  has  the  proper  form 
for  t  — >0  is 


p  =  (x  +  UQt  -  q)/at 


where 


q  =  i/2 


/2  (xo  +  abt) 

[-(‘o  *  4  - 


This  mathematical  solution  is  physical¬ 
ly  valid  for  0  <  t  <  (abxQ)  ^  and 

0<x  +  ut<x  . 

o  o 

We  generate  the  negative  lobe  by  fold¬ 
ing: 


p(-x-uQt,  t)  =  -  p(x+uQt,  t).  (20) 


For  t>(abx  )  ,  the  mathematical  solu- 
o 

tion  is  multivalued,  with  Eq.  (19) 
joining  at  a  branch  point  onto  the  cor¬ 
responding  root  with  positive  radical. 
The  physical  solution,  however,  follows 
Eq.  (19)  out  to  the  point  at  which  the 
area  under  the  positive  lobe  equals 
that  of  the  initial  condition,  and  then 


f  XS  /»  Xp 

/ p(x,t) d x  =  /  p(x,o)  dx 
*^-u  t  *'o 


From  Eqs.  (14),  (18),  and  (19)  this 
gives 


L  =  _  w  +  _l\ 

it  2  2  \  o  abt / 

-  f  abt  [?  ("o  +  St)  -S?]^  122  > 
Wx  +_l\3 

12  \xo  abt/  ’ 


where  x’  =  x  +  u  t. 

o 


Even  though  x'=xq  satisfies  Eq .  (21) 

for  all  t>0,  it  does  not  give  the  shock 
location.  This  root  holds  for  the 
steepening  phase  of  the  physical  solu¬ 
tion,  and  after  shock  formation  it 
applies  only  to  the  double  valued  math¬ 
ematical  solution.  Shock  formation  oc¬ 
curs  when  the  radical  in  Eq.  (19)  first 

vanishes  on  the  interval  0<x' <x  .  This 

o 

occurs  at  (x',t)  =  (x  ,  (abx  )  S.  Eq. 

o  o 

(22)  yields  a  quartic  equation  for  x', 
that  may  be  reduced  to  a  quadratic  by 

removal  of  the  double  root  x'”x  .  Cf 

o 

the  two  remaining  roots,  the  one  cor¬ 
responding  to  a  rightward  progressing 
shock  is 
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saving  computer  time  and  avoiding  the 
risk  of  a  zero  divide. 


-uot 


+ 


9abt 


(24) 


In  summary,  the  parabolic  lobe  analytic 

solution  consists  of:  P =0  for  x>x  ; 

s 

Eqs.  (18)  and  (19)  for  -u  t<x<x  ;  and 

o  s 

Eq.  (20)  for  x<-u  t. 


4.1  The  First-Order  Scheme 

A  representation  of  Eq.  (la),  which  is 
first-order  and  which  reduces  to  the 
familiar  upwind  scheme  in  the  case  of 
constant  flow  velocity,  is  as  follows: 


f . 
1 


+  1/2 


+  f. 


1/2, 


(27) 


The  second  analytic  solution  to  be  used 
in  numerical  tests  is  that  of  a  square 
wave  moving  at  a  constant  velocity  u  : 


1  ,  0<x  +  uQt<cx0 
0  otherwise. 


(25) 


This  simple  case  will  serve  to  illus¬ 
trate  effects  of  diffusion  and  disper¬ 
sion  in  the  numerical  algorithms. 

4.0  Numerical  Techniques 

We  present  two  separate  upwind  differ¬ 
ence  schemes  and  combine  them  in  a  hy¬ 
brid  scheme  by  the  method  of  flux  cor¬ 
rection  (Zalesak,  1979).  For  simplicity 
we  assume  a  constant  grid  interval  fix. 
We  take  and  4^  to  be  values  at  inte¬ 
ger  grid  points,  and  define  half  inte¬ 
ger  grid  point  variables  f^  +  an4 

Wi  +  1/2*  ®ere  f  l®  a  numerical  flux  to 

be  specified  below  and  w  is  the  sign  of 
the  characteristic  speed  (whose  only 
purpose  is  to  define  the  upwind  direc¬ 
tion).  In  practice  we  take 

wi  +  1/2  =  (*i+l  '  ♦i)*(p1+l  ‘  pi}-  (26) 

Note  we  have  replaced  the  division  in  a 
second  order  representation  u  in  Eq. 
(3)  by  a  multiplication.  This  results 
in  w  being  of  the  proper  sign,  while 


where  the  superscript  indicates  the 
time  level  (t  =  n6t)  and 


fi  +  1/2 


V2  58  °  (28) 


4^  6t/6x  ,  W.  + 

4i+1  <St/6x  ,  Wi  +  1/2  <  °  * 


The  distinction  between  this  scheme  and 
ordinary  upwind  lies  in  the  explicit 
prescription  Eqs.  (26)-(28)  for  the 
upstream  direction  in  terms  of  point 
values  of  <(>  and  p  as  opposed  to  the 
characteristic  speed  u  or  the  flow 
speed  v. 

Quantities  without  superscripts  are 
evaluated  at  time  level  n.  The  above 
scheme  preserves  monotonicity,  as  shown 
by  consideration  of  the  following  four 
cases : 


Case  1.  w..  _  1^2  ?0,  w..  +  ^  >0. 


Eqs.  (27)  and  (28)  give 


P^1  -  p”  *  -  6t/6x  (<(>..  -  4>^_-j)  (29) 


-  6t/6x  u_(p.  -p._.|)  (30) 


where 


and 


U-  -  {$i  -  4>i  -J  )/(pi  -  P^-J ) 

=  u(p.  ). 


(31) 


u+  =  ^i  +  1  ‘  *i)/(^i+l  "  •' 0 

=  u(p+). 


(36) 


The  last  identity  results  from  the  mean 
value  theorem,  with  p.  lying  In  the  in¬ 
terval  (p^,  Then 


with  P  +  lying  on  the  interval  (p  p 
P^).  Monotonicity  is  preserved  in  this 


case  for  e+<l . 


p”+1  =  Pj  (1  -  e_)  +  P^-j  e_  »  (32) 


Case  3.  w.  _  y2  <  0,  w.  +  ]/2  >  0 


where 


Eqs.  (27)  and  (28)  yield 


e_  =  | U_  >t/6x|  . 


(33) 


n+1 


(37) 


Even  though  u-  is  positive  in  this 
case,  we  retain  the  absolute  value  in 
Eq .  (33)  for  convenience  below.  Eq . 
(32)  demonstrates  that  monotonicity  is 
preserved  for  c^l.  Specifically,  the 
right-hand  side  of  Eq .  (32)  consists  of 
a  weighted  mean  of  neighboring  p  val¬ 
ues,  with  each  weight  being  non¬ 
negative  for  e<l. 


Case 


2.  w._1/2  <  0,  wi+1/2  < 


0. 


As  above,  we  find 


n+1 


i  =  P-j  (]  "  e+)  +  e-j+i  £+  ’ 


This  reflects  that  u=0  within  one  cell 
of  point  i. 


Case  4-  wi  _  y2  >0,  wi  +  1/2  <  0 


The  straightforward  result  from  Eqs. 
(27)  and  (28)  gives 


n+1 


Pi  =  -6t/6x  (<f>i+1  -  t^j).  (38) 


A  small  amount  of  manipulation  after 
adding  and  subtracting  <J>  within  the 

parentheses  yields 


where 


e+  =  |u+  6 1 / 6  X |  , 


(35) 


n+1  _  ,,  x  , 

p.  -  p.  (1-t.  -  e+)  +  „i+1  t+ 

>1-1 


(39) 
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This  case  provides  the  most  stringent 
condition  on  monotonicity  (and  thus 
stability  as  well): 


| u  6t/5x|  <  1/2  .  (40) 

4.2  The  Higher-Order  Scheme 

A  one-sided  representation  of  Eq.  (la) 
that  has  second-order  spatial  accuracy 
is 


where 


Si  +  1/2 


fi  -  1/2  ’  Wi  +  1/2  >0 
fi  +  3/2  5  wi  +  1/2  <0‘ 


(45) 


The  generalization  of  Eq .  (41)  is  then 


p"+1  -  pj  =  -  6t/6x  (3/2  <(,.  -  2<j>i._1 


+  l/24>._2). 


(41) 


This  representation  has  higher-order 
accuracy  than  the  first-order  scheme 
Eqs.  (26)-(28),  but  is  of  little  in¬ 
terest  by  itself  because  it  does  not 
preserve  monotonicity.  Under  the  as¬ 
sumption  that  the  characteristic  speed 
is  locally  constant,  the  stability  con¬ 
dition  for  Eq.  (41)  at  high  wavenumbers 
is 


"Fi+l/2  +  Fi  -  1/2 


(46) 


In  practice,  F  is  constructed  according 
to  Eqs.  (44)  and  (45),  but  Eq .  (46)  is 
replaced  by  the  hybrid  scheme  to  fol¬ 
low. 


4.3  The  Hybrid  Scheme 

Following  Zalesak  (1979)  we  construct 
the  following  provisional  quantities: 


"  f i  +  1/2 


+  f. 


1/2 


(47) 


|  u  6t/6x |  s  1/2  ,  (42) 

which  coincides  with  Eq .  (40).  Note 
that  the  first-order  analog  of  Eq.  (41) 
is 


p"+1  -  p"  =  -  <5t/6x  (<pi  -  )  .  (43) 


By  comparison  of  Eqs.  (43)  and  (41)  we 
construct  numerical  fluxes  for  the 
second-order  scheme  as  follows: 

Fi  +  1/2  =  3/2  fi  +  1/2  '  1/2  Si  +  1/2  5 

(44) 


6fi  +  1/2  =  Fi  +  1/2 


fi  +  1/2  • 


(48) 


We  now  have  a  first-order  update  C*  in 
which  raonotonicity  has  been  preserved. 
Next,  we  perform  a  "flux  correction" 
before  completing  the  update.  This  is 
simply  a  filtering  operation  that  for¬ 
bids  higher-order  enhancement  of  extre¬ 
ma.  For  simplicity  we  use  the  original 
form  of  Boris  and  Book  (1973): 


6fi  +  1/2  -  s 

[l  <5fi  +  1  /2  i  ’ 

Si  +  1/2  (pi 


i  +  1/2 


max 


^0,  min 


c  (  .  *  ) 

5i  +  1/2  u  i+2  ’■'i  +  l  ’ 


(49) 
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where 


Si  +  1/2 


3fi  +  1/2^° 
5fi  +  1/2  <0 


(50) 


Innovations  to  the  above,  including 
extension  to  multidimensions  without 
timestep  splitting,  have  been  given  by 
Zalesak  (1979). 


where  the  flow  speed,  rather  than  the 
characteristic  speed,  was  used  to  de¬ 
termine  the  upwind  direction.  The  third 
case  is  the  passive  advection  of  square 
wave  profile.  This  example  illustrates 
most  clearly  the  gains  achieved  through 
higher  order  flux  correction.  All  exam¬ 
ples  use  a  mesh  of  202  points  spaced  at 
equal  intervals.  The  two  endpoints  are 
reset  after  each  calculation  using  per¬ 
iodic  boundary  conditions.  The  timestep 
is  constant  for  each  case,  resulting  in 
Courant  numbers  u  St/Sx  as  given 
below. 


The  hybrid  scheme  is  completed  as  fol¬ 
lows  : 


»r  *  -  5fi  *  1/2  ♦  sfi .  1/2  -<5i> 

One  advantage  of  the  above  hybrid 
scheme  over  schemes  involving 
piecewlse-analytic  solutions  is  its 
programming  simplicity.  The  subroutine 
that  executes  the  entire  hybrid  scheme, 
Eqs.  (27-51),  consists  of  29  executable 
Fortran  statements,  of  which  six  are 
calls  to  an  endpoint  setting  routine  of 
length  six  statements  (see  the  Appen¬ 
dix).  A  second  advantage  is  that  the 
scheme  itself  is  independent  of  the 
specific  relation  between  <J>  and  p 
(i.e..  Burgers'  equation  or  whatever). 

5.0  Numerical  Results 

In  this  section  we  present  calculations 
for  three  examples  in  which  exact  re¬ 
sults  are  known.  Each  is  based  on  re¬ 
sults  given  in  Section  3.  The  first  is 
the  development  of  an  N-wave  shock  in 
which  positive  and  negative  lobes  pro¬ 
duce  an  antisymmetric  velocity  distri¬ 
bution  (i.e.,  uQ  "  0  in  Eq.  (13)).  Sec¬ 
ond  is  the  development  of  an  N  wave  in 
drifting  frame  (uq*0).  The  motivation 

for  this  case  is  to  demonstrate  that 
expansion  shocks  do  not  result  where  u 
■  0,  p  #0  and  ux>0.  We  noted  that  they 

occur  under  these  circumstances 


Example  1. 

Initial  conditions  are  taken  from  Eqs. 
(14)  and  (15).  Time  advancement  is 
achieved  by  defining 

*1  =  a/2  of  ~  u0ci  (52) 


as  in  Eq .  (13)  (with  uq  =  0  in  this 

case)  prior  to  executing  the  hybrid 
scheme  of  Section  4.  The  timestep  is 
such  that  the  maximum  Courant  number  is 

ab  xQ2  6 1 / 46x  =  1/4.  (53) 


Figure  1  shows  the  Initial  condition 
followed  by  results  from  the  first- 
order  scheme  and  results  from  the  hy¬ 
brid  scheme.  Parameters  are  such  that 
the  shock  forms  at  step  41.  After  this 
point  in  time,  the  second-order  correc¬ 
tion  shows  increasing  improvement  over 
the  purely  first-order  result.  Note 
that  the  second-order  method  at  step 
601  shows  clear  superiority  in  locating 
the  shock  properly  and  maintaining  its 
sharpness . 

Example  2. 

Initial  conditions  are  the  same  as  in 

example  1,  but  a  drift  term  u  ^  0  is 

o 

present  in  Eq.  (52).  The  timestep  is 
the  same  as  in  example  1,  but  the 


did 


6.0  Conclusions 


presence  of  the  drift  term  raises  the 
Courant  number  from  0.25  to  0.3375.  The 
value  of  u  is  such  that  u  “  0,  where 

is  equal  to  0.35  of  its  initial  maximum 
value.  As  shown  in  Figure  2,  the  pro¬ 
files  remain  smooth  at  this  value,  dem¬ 
onstrating  stability  against  nonphys¬ 
ical  expansion  shocks.  Tests  run  with 
higher  values  of  uq  have  also  confirmed 

this  point.  In  tests  using  the  local 
flow  speed  V  =  l/2ap  -  UQ  rather  than  u 

to  determine  the  upwind  direction,  ex¬ 
pansion  shocks  did  occur  where  v  =  0 
and  v^  >  0.  Comparison  of  the  first- 

order  portions  of  Figures  1  and  2  il¬ 
lustrates  the  first-order  scheme's  in¬ 
creasing  dissipation  with  increasing 
Courant  number.  The  second-order  re¬ 
sults,  however,  are  much  less  subject 
to  this  type  of  error.  The  steplike 
structure  at  step  601  is  simply  a  re¬ 
sult  of  the  periodic  boundary  condi¬ 
tions  used  here. 


Example  3. 

The  square  wave  profile  shown  in  Figure 
3  is  generated  from  Eq.  (25).  This  test 
consists  of  linear  advection  with  a 
courant  number  of  1/4.  At  the  beginning 
of  each  step  the  flux  is  defined  from 
Eq.  (52)  with  a  •  0.  As  the  calculation 
proceeds,  the  first  order  results  of 
Figure  3  show  the  initially  sharp 
structure  being  progressively  smoothed 
out.  The  second-order  results,  however, 
retain  a  much  greater  degree  of  fidel¬ 
ity  to  the  analytic  solution. 

The  results  of  Zalesak  (1981)  suggest 
that  the  residual  dissipation  present 
in  the  second-order  results  could  be 
further  reduced  by  employing  higher- 
order  upwind  differences  in  the  correc¬ 
tion  stage  (44)-(51)  of  the  hybrid 
scheme . 


We  have  presented  a  simple  first-order 
scheme  for  solving  Eq.  (la)  which  main¬ 
tains  montonicity  and  avoids  artificial 
expansion  shock  generation  where  the 
characteristic  speed  reverses  sign. 
This  is  done  without  the  use  of  arti¬ 
ficial  viscosity  or  Riemann  solutions. 
We  have  also  shown  that  the  use  of 
second-order  upwind  differences  and  the 
method  of  flux  correction  greatly  re¬ 
duce  the  dissipative  effects  of  the 
first-order  method  while  preserving 
monotonicity.  The  Rankine-Hugoniot  re¬ 
lation  for  shock  propagation  speed  is 
implicit  in  the  method  by  virtue  of 
flux  conservation  at  the  shock  front. 
Comparisons  of  numerical  and  analytic 
results  for  a  developing  N-wave  shock 
and  for  linear  advection  of  a  square 
wave  confirm  the  accuracy  of  the 
method. 
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INITIAL  CONDITION 


STEP  I 


Figure  2.  Drifting  N-wave  demonstrating  stabl I 
Ity  against  expansion  shocks.  The  characteris¬ 
tic  speed  vanishes  where  the  amplitude  is  35 
percent  of  Initial  maximum.  For  legend  see 
Figure  1.  The  step-1  Ike  structure  at  step  601 
Is  a  result  of  periodic  boundary  conditions. 


INITIAL  CONDITION 


STEP  1 


Figure  3.  Linear  advectlon  of  a  square  wave. 
Solid  line:  exact  solution;  dashed  line: 
numerical  solution.  First-order  results  appear 
In  b-d;  second  In  e-g. 


Appendix:  Fortran  Statement  Listing  of 
the  Upwind  Flux  Method  Subroutine 


t 


i 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 

c 


c 

c 


c 


SUBROUTINE  FLUX(N.  PHI.  F.S.U.  R> 

REAL  F'HI(N),  F(N>.  S(N).  W<N>.  R  <  N  ) 

VERSION  OF  MAY  19e3 

INTEGRATE  FORWARD  ONE  TINESTEP:  DR/DT  =  -D  PHI/DX 

1ST  ORDER  UPWIND  ...  FLUX  CORRECTED  TO  2ND  ORDER. 

PHI  HAS  BEEN  NORMALIZED  BY  DT/DX  IN  THE  CALLING  PROGRAM. 

SHORT  WAVELENGTH  STABILITY  CRITERION  IS  U  ,5. 

WHERE  U  =  D  PHI/DR  IS  THE  CHARACTERISTIC  SPEED. 

R  AND  PHI  RESIDE  AT  INTEGER  GRID  POINTS  AND  MUST  BE  PREDEFINED. 
FLUX  VARIABLES  FiS.W  ARE  SCRATCH  ARRAYS  RESIDING  AT  HALF  INTEGER 
GRID  POINTS.  THE  INDEXING  CONVENTION  IS 

F<I)  RESIDES  AT  I-.5  ...  SAME  FOR  S  AND  W. 

Nf  =  N-l 

DO  2  I=2.NM 

W(I)  =  ( F'HI  ( I  )  -  PHI  ( 1-1  >  >*<R(  I )  -  R  <  I  - 1  )  ) 

F ( I )  =  PHI(I-l) 

IF(U( I )  . LT .  0. )  F ( I )  =  PHI < I ) 

2  CONTINUE 

CALL  PBCSET ( N  .  F> 

DO  3  I =2  *  NM 
1ST  ORDER  UPDATE 

R  ( I  )  =  R  <  I  )  -  F  (  I  + 1  )  +  F  (  I  ) 

PREPARE  2ND  ORDER  CORRECTION 
FI  =  F(I-l) 

IF  <  U  < I )  .LT .  0. )  FI  =  F( 1  +  1  ) 

3  S  < I >  =  .5*<F< I >  -  FI) 

CALL  PBCSET(N.S) 

CALL  PBCSET  <  N . R ) 
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FLUX  LIMITER 

DO  A  1=2. NM 
F ( I )  *  SIGN< 1 . .  S  < I >  ) 
U  ( I  )  =  R  <  I  )  -  R  (  I  -  I  ) 
CALL  PBCSET(N.W) 


DO  5  I-2.NM 

S(I)  =  F< I >*AMAX1 (0. .AMIN1 < ABS(S< I ) ) tF< I )*U< 1  +  1 ) .F< I )»U<  1-1  )  )  ) 
CALL  PBCSETCN.  S> 


FINAL  CORRECTION 
DO  6  1=2. NM 

6  R< I )  =  R< I )  -  S< 1  +  1 )  +  SI  I ) 
CALL  FBCSET(N.  R) 

RETURN 

END 


SUBROUTINE  PBCSET(N.  A) 
REAL  A ( N ) 

SET  PERIODIC  BOUNDARIES 
A (  1  )  =  A  <  N- 1 ) 

A  (  N  (  =  A  (  2  ) 

RETURN 
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